arXiv: 1508.03189vl [cond-mat.mes-hall] 13 Aug 2015 


Merging of the Dirac points in electronic artificial graphene 


J. Feilhauer, 1,2, 0 W. Apel, 1 and L. Schweitzer 1 

1 Physikalisch-Technische Bundesanstalt (PTB), Braunschweig, Germany 
2 Institute of Electrical Engineering, Slovak Academy of Sciences, 841 04 Bratislava, Slovakia 

(Dated: August 14, 2015) 

Theory predicts that graphene under uniaxial compressive strain in an armchair direction should undergo a 
topological phase transition from a semimetal into an insulator. Due to the change of the hopping integrals 
under compression, both Dirac points shift away from the corners of the Brillouin zone towards each other. For 
sufficiently large strain, the Dirac points merge and an energy gap appears. However, such a topological phase 
transition has not yet been observed in normal graphene (due to its large stiffness) neither in any other electronic 
system. We show numerically and analytically that such a merging of the Dirac points can be observed in 
electronic artificial graphene created from a two-dimensional electron gas by application of a triangular lattice 
of repulsive antidots. Here, the effect of strain is modeled by tuning the distance between the repulsive potentials 
along the armchair direction. Our results show that the merging of the Dirac points should be observable in a 
recent experiment with molecular graphene. 

PACS numbers: 73.20.At, 73.22.Pr 


I. INTRODUCTION 

In the past decade, graphene^ became one of the most stud¬ 
ied materials due to its unusual properties that are attractive 
from the perspective of basic physics as well as from that 
of a future technology. In graphene, the two-dimensional 
hexagonal lattice of carbon atoms provides for a gapless elec¬ 
tronic bandstructure which consist of two bands touching each 
other at two Dirac points located at the independent corners 
of the hexagonal Brillouin zone. The Fermi energy of natu¬ 
ral graphene lies in the vicinity of the Dirac points, where the 
bands have linear dispersion and conical shape. Therefore at 
low energies, the bandstructure of graphene can be modeled 
by two Dirac cones and electrons behave as massless Dirac 
particles with a valley isospin. 

After the discovery of graphene, there has been a significant 
effort 277 — to create artificial systems with hexagonal symme¬ 
try that mimic the unique properties of Dirac quasiparticles. 
The successful experimental realization of Dirac quasiparti¬ 
cles was reported in systems of cold atoms in hexagonal op¬ 
tical lattices^, in photonic crystals created by a hexagonal ar¬ 
ray of optical waveguides^ and in two-dimensional electron 
gases with hexagonal supcrlattices^— The main advantage of 
these artificial systems is that one can control and tune the sys¬ 
tem parameters (i.e., the lattice constant and hopping param¬ 
eters) independently. That makes it possible to study features 
that are not observable in normal graphene, e.g., the extremely 
large pseudomagnetic fields induced by strain, Kekule lattices, 
edge states, etc. 

One of the interesting possibilities in artificial graphene is 
to induce lattice anisotropy by tuning the hopping parameters 
between the lattice sites. The theory based on the nearest- 
neighbor tight-binding approximation predicts^— that by 
changing the hopping parameters, both Dirac cones are shifted 
away from the corners of the Brillouin zone and become el¬ 
liptical instead of circular. With increasing anisotropy, the 
Dirac cones merge and an energy gap is created. This causes a 
topological phase transition-^ 7 where the system turns from a 
semimetal into an insulator. At the critical point of this merg¬ 


ing, the energy bands touch only in a single point with an 
unusual dispersion relation in its vicinity.— In one direction in 
k-space, the dispersion is linear as in the standard Dirac cone 
but in the perpendicular direction, the dispersion is quadratic 
as in the case of massive electrons. 

In normal graphene, lattice anisotropy could be induced by 
application of uniaxial strain that changes the distance be¬ 
tween the lattice sites. But the values of strain that would lead 
to the topological phase transition are unrealistically large ow¬ 
ing to the large graphene stiffness^ Thus, merging of the 
Dirac points was demonstrated experimentally in systems of 
cold atom ‘4-12 and microwave photonic crystals^ where the 
distance between lattice sites can be varied much more than 
in normal graphene. In spite of some efforts^ a merging of 
the Dirac points has not yet been confirmed in any electronic 
system. 

Here, we show analytically and numerically that a merging 
of Dirac points appears also in electronic artificial graphene 
under uniaxial compressive ’strain’. We study artificial 
graphene created in a two-dimensional electron gas modulated 
by a triangular lattice of repulsive potentials, e.g., repulsive 
molecules on a copper surface^ or a muffin-tin potential of 
antidots in a semiconductor heterostructure^— In such sys¬ 
tems, the electrons are repelled into the space between the 
antidots which leads to hexagonal symmetry. The lowest en¬ 
ergy bands of artificial graphene are similar to the two-level 
tight-binding bandstructure of natural graphene although the 
electrons are not bound to any attractive potential. The lat¬ 
tice constant of artificial graphene can be orders of magnitude 
larger than in natural graphene and, moreover, the position of 
the artificial lattice sites can be easily tuned. Strain is simu¬ 
lated by tuning the distance between the repulsive triangular 
potentials along the armchair direction. When the antidot lat¬ 
tice is stretched, the Dirac points are shifted away from the 
corners of the Brillouin zone into its interior, but they always 
exist. On the other hand, with increasing compressive strain, 
both Dirac points are moving along the edge of the Brillouin 
zone towards each other until they merge with a dispersion as 
described above. For larger compression, a bulk gap is ere- 
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ated. 

We believe that our theoretical results could be verified ex¬ 
perimentally in the recently created molecular graphene,^ 
where a two-dimensional electron gas on the Cu surface 
is modulated by a triangular antidot lattice of repulsive 
molecules (CO or coronene). Here, the presence of a Dirac 
energy spectrum is proven by measurements of the density of 
states (DOS) which is similar to typical graphene DOS with 
the zero value at the energy corresponding to the Dirac point. 
The authors of performed an experiment where the artifi¬ 
cial graphene was elongated in the armchair direction by about 
30%. They observed no gap opening in the DOS and that is 
consistent with our numerical results for the case of streched 
artificial graphene. Nevertheless, our numerical calculations 
of the DOS show that the molecular graphene with coronene 7 
is a very promising candidate to observe the merging of the 
Dirac points and related gap opening in the DOS. This should 
occur, however, for a lattice compression of about 25%, and 
that is technologically feasible. 

This paper is organized as follows. In the second section, 
we review the results of the analytical tight-binding calcu¬ 
lations that describe the merging of the Dirac points in the 
hexagonal lattice with the hopping matrix anisotropy. In the 
third section, we define mathematically the model of artificial 
graphene and discuss the effect of uniaxial strain. Then we 
describe the technique that we used to calculate the numerical 
results. The fourth section contains our numerical results of 
the bandstructure for artificial graphene under uniaxial strain. 
The calculations show that the merging of the Dirac points 
occurs in the artificial graphene under compressive strain in 
the armchair direction. The fifth section shows the analyt¬ 
ical calculation of the low-energy bandstructure of strained 
artificial graphene in the nearly-free electron approximation, 
which is valid for weak potentials. The analytical formulas 
confirm the merging of the Dirac points in the compressed 
artificial graphene and are in a good agreement with the nu¬ 
merical results. In the sixth section, we discuss the realiz¬ 
ability of our calculations in real experimental situations and 
we conclude that the merging of the Dirac points should be 
observable in the recent experiment with molecular graphene 
with coronene . 7 


II. MERGING OF THE DIRAC POINTS: TIGHT-BINDING 
MODEL DESCRIPTION 

In this section, we review the results of the nearest-neighbor 
tight-binding calculations^— concerning the bandstructure 
of a hexagonal lattice with anisotropy in the hopping matrix 
elements. We examine the positions of the Dirac points in the 
Brillouin zone. 

Normal graphene consists of a triangular Bravais lattice 
with a pair of carbon atoms located in its primitive cell. The 
lattice vectors are 

u _ £„ V3£_ u £_ , 

t>i — ^ e x ^ e y> * 32 — 2 ex + 2 e ^’ 1 

where L is a lattice constant and Lj\J 3 is the distance be- 



FIG. 1: (a) Model of the hexagonal lattice of graphene. Each lat¬ 
tice site has three nearest-neighbors with the locations given by the 
vectors 5i (see Eq[3]t and corresponding hopping parameters ti. (b) 
The repeated Brillouin zone scheme of graphene. The Brillouin zone 
contains a symmetry point T and a pair of independent corners K 
(green empty circles) and K’ (green full circles) which are period¬ 
ically repeated around the reciprocal lattice points (black crosses). 
When t\ ■ t:i — £ and t 2 = t', the Dirac points are located on a 
line CC'. For 0 < t' /f < 1, the Dirac points are located inside the 
Brillouin zone (red lines) and for 1 < t'/t < 2 the Dirac points lie 
at the edge (blues lines). For t'/t = 2 the Dirac points merge at S 
(black dots). 


tween neighboring carbon atoms. The corresponding recip¬ 
rocal lattice is triangular as well with a hexagonal Brillouin 
zone with the side length 47 r/( 3 L). In the nearest-neighbor 
tight-binding approximation, each carbon atom possesses one 
localized electronic state that slightly overlaps with its three 
nearest neighbors. These overlaps are in general character¬ 
ized by three hopping parameters t\, t- 2 , (see Fig. [Qa)). 
In the isotropic case, all hopping parameters are equal. If 
we take the hopping anisotropy in armchair direction, we get 
ft = f 3 / t- 2 - Then, the bandstructure can be written in the 
form 

e(k) = ±| te tkSl + t'/3e ikS2 + te lkS3 \, ( 2 ) 
where we have defined fi = £3 = t, t 2 = t' and 

{<5i, <52, £ 3 } = — {—bi — 2b2, b 2 — bi, 2bi + b 2 } , (3) 

are the vectors connecting an atom with its neighbors. 

The bandstructure © consists of two bands that touch at 
zero energy at the Dirac points. We obtain their position in 
the reciprocal space by solving e(k) = 0 which implies the 
conditions 

c°s (^Y Lk y^J =±l, cos (^Lk^=^ f —. (4) 
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The conditions (|4|) give two independent Dirac points per Bril- 
louin zone which are periodically repeated through the whole 
reciprocal space. To study the evolution of the Dirac points 
with increasing hopping anisotropy, it is more convenient to 
use the repeated Brillouin zone scheme and to focus on the 
Dirac points Di and D 2 located at 
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Di _ J.D; 


fc? 1 = 2 
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For the case of hopping anisotropy in the armchair direction, 
the position of the Dirac points in reciprocal space is restricted 
to the k x axis. Moreover, the positions of Di and D 2 are sym¬ 
metric around the point S at k x = 2ix/L. Therefore, in the 
following text, we discuss mostly only the position of Di. 

The position and existence of Dirac points is given by the 
value of the anisotropy parameter t'/t. Firstly, if 0 < t'/t < 
1, the Dirac point Di lies inside the Brillouin zone (red lines in 
Fig.HJb)) and the lowest value of fcP 1 is limited by the value 
7 t/L for t'/t = 0 (point C). For t'/t = 1 (isotropic lattice), 
the Dirac point lies exactly at the corner of the Brillouin zone 
K with fcP 1 = 47 t/3 L. Secondly, if 1 < t'/t < 2 the Dirac 
point is located at the edge of the Brillouin zone (blue lines in 
Fig.HJb)) and for the critical point t'/t = 2, both Di and D 2 
merge at the point S. Thirdly, for 2 < t'/t, Dirac points no 
longer exist and an energy gap opens. 

The hopping anisotropy discussed above could be achieved 
in normal graphene by applying an uniaxial strain in the arm¬ 
chair direction. When the graphene sample is stretched, the 
interatomic distance | ^ 2 1 increases. Then, the overlap of elec¬ 
tron wavefunctions is suppressed and the hopping parameter 
(2 decreases. This corresponds to the situation where t'/t < 1. 
On the other hand, when the sample is compressed, the inter¬ 
atomic distance |<5 2 | decreases and (2 increases. Then t'/t > 1 
and it increases with increasing compression. In principle, it 
could be possible to observe a merging of the Dirac points (for 
t'/t = 2) in graphene, but, unfortunately, the mechanical stiff¬ 
ness of graphene is very larged and the sufficient compression 
to obtain t'/t = 2 is unreachable. 

To observe a merging of the Dirac points, one has to ar¬ 
range a system with hexagonal symmetry where it is possible 
to tune the hopping parameters in a wider interval than in nor¬ 
mal graphene. This was already shown in a lattice of cold 
atoms2ii£ or a photonic crystal^ but not yet for any electronic 
system. Our aim is to study a merging of the Dirac points in 
electronic artificial graphene described in the next section. 


FIG. 2: (a) Triangular lattice of antidot potential spanned between the 
vectors ai and a 2 . The antidot radius is p/L = 0.2 and the lattice 
is compressed in the y (armchair) direction with a = 1.25. (b) The 
corresponding reciprocal lattice (crosses) is triangular as well and it 
is strained in the horizontal direction. In general, the first Brillouin 
zone is an irregular hexagon and we use the repeated Brillouin zone 
scheme (dashed lines). Because the antidot lattice is compressed, the 
Dirac points Di and D 2 are shifted away from the corners Ki and 
of the Brillouin zone towards the point S. 

is to introduce a hexagonal lattice of potential wells, which 
was performed lithographically in a GaAs heterostructure^ 
In the present paper, we study an artificial graphene created 
from a two-dimensional electron gas by application of a repul¬ 
sive triangular potential of antidots V(r) (see Fig. [2] a)). This 
kind of artificial graphene was extensively studied^ - - 23 and 
it captures the main features of the experiments with molecu¬ 
lar graphene, which is a two-dimensional electron gas on the 
copper surface modulated by the repulsive molecules of CQ& 
and coronene. 7 

It can be shown& that the wavefunctions of the electrons 
from the two lowest energy bands are localized around the 
centers of the triangles formed by the antidots and thus form 
a hexagonal lattice. Consequently, the two lowest energy 
bands of artificial graphene are very similar to those of nor¬ 
mal graphene. Namely, there is a pair of Dirac points with 
linear dispersion in their vicinity. 

The effect of strain discussed in the previous section can 
be easily achieved in this artificial graphene by prolonging or 
reducing the dimension of the triangular antidot lattice in the 
armchair direction. This consequently changes the distances 
in the underlying hexagonal electron lattice. 

A. Trigonal antidot potential under uniaxial strain 


III. ARTIFICIAL GRAPHENE UNDER UNIAXIAL 
STRAIN 

There have been several theoretical and experimental at¬ 
tempts to create an electronic artificial graphene. Most of 
them are based on a nanopatterning of a two-dimensional elec¬ 
tron gas by an external periodic potential. One possibility 


The antidot potential V (r) consists of repulsive circles with 
radius p ordered on a triangular lattice with lattice constant L 
and reads 

U(r)=V^0(p-|r-R|), (8) 

R 

where V is the potential strength, 0 is the Heaviside step func¬ 
tion and R = roai + n &2 are the lattice sites spanned by the 
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the Bloch form 
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(9) 


V'k(r)=^c^e i ( k+G ) r , (15) 

G 


where the parameter a represents an uniaxial ’strain' in the y 
(armchair) direction. For a < 1, the lattice is prolonged in 
the y direction and for a > 1 it is compressed. In the case 
a = 1, we obtain the unstrained triangular lattice vectors ©. 
The corresponding reciprocal lattice to I© is G = fca* + /a?j, 
where 


a 
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( 10 ) 


The Brillouin zone of such a lattice is hexagonal (see 
Fig. ©b)), but in contrast to the calculations in section II, it 
is not a regular hexagon unless a = 1. Its side lengths are 
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and the positions of Ki and K/ points are 




47T 

T 



( 12 ) 


We would like to stress that the uniaxial constant strain ap¬ 
plied to the antidot lattice does not act as a constant strain on 
the underlying hexagonal lattice of electrons. For the strain 
a = v/3 « 1.73, the triangular antidot lattice becomes square 
and the corresponding ’electron’ lattice changes from hexag¬ 
onal to square. A constant strain on the hexagonal lattice 
would only lead to the deformation of the regular hexagons 
but would not change the connectivity of the lattice. For 
a = v/3 also the Brillouin zone is no longer hexagonal but 
square and we get |Ki — K'J = 0. For a > v/3, the antidot 
lattice becomes the same as the lattice for a < x/3 but it is 
rotated by 90° with the rescaled lattice constant L' = a 
and strain parameter a' = 3/a. For example, for a = 3, 
the antidot lattice again corresponds to the unstrained artifi¬ 
cial graphene with L' = L/v/3 and zig-zag orientation in the 
y direction. 

The previous discussion shows, that it is sufficient to restrict 
the values of a to the interval [0, -\/3] - 


where k is restricted to the first Brillouin zone. The 
Schrodinger equation (IT3l) then becomes 


E 


'h 2 { k + G) 2 
2 TO* 


Sg,,g' + Vq' 



L.Q/ 


o, (16) 


where the Fourier coefficients of the lattice potential are 

yQ = ^ (17) 

where Ji is the Bessel function of the first kind and S = 
\/3L 2 / (2a) is the area of the unit cell. Eq. (ITTb shows that 
due to the rotational symmetry of an antidot, the Fourier coef¬ 
ficients of V (r) depend only on the length of the wavevector 
Q. Eq. dl~6b shows that the periodic lattice potential mixes 
only those plane waves with wavevectors that are shifted by 
reciprocal lattice vectors G. Solving (fl6l > numerically we get 
the complete band structure e„(k). 

In general, the number of terms in the sum ( fl~5T > is infinite 
(and the number of the corresponding equations (IT6l > as well) 
but in practice it is sufficient to assume a finite number of re¬ 
ciprocal lattice vectors G. Because we are interested in the 
lowest energy bands e n (k), the most important terms in (fl5l i 
are the plane waves with the smallest wavevectors. Therefore, 
we restrict the summation in ( fiTl i to the wavevector F = 0 
plus a few hexagonal shells of reciprocal lattice points cen¬ 
tered around T. The calculated energy bands saturate quickly 
with increasing number of hexagonal shells N and a sufficient 
value of N depends on the strength of potential. In our calcu¬ 
lations, we use mostly N = 20. 

It can be easily shown that the electron eigenenergies in 
(fl~6t scale as h 2 /(to* L 2 ) and therefore, we introduce the en¬ 
ergy unit 


h 2 ( 47 t \ 2 
2m* \3 LJ 


(18) 


which corresponds to the free electron energy at the K point. 
Then, we use the dimensionless parameters E/Eq, V/Eq and 
p/L. 


B. Numerical model 

The electron wavefunction in artificial graphene is de¬ 
scribed by the Schrodinger equation 

Hip(T) = eip(r), (13) 

with the Hamiltonian 

H '-2^ A + F(r) ' (14) 

where the lattice potential V (r) given by © is periodic with 
the lattice vectors R. Then the electron wavefunctions are of 


IV. NUMERICAL RESULTS 

In this section, we present the results of our numerical cal¬ 
culations describing the bandstructure of artificial graphene 
under uniaxial strain. Firstly, we show the results for the un¬ 
strained case a = 1. 


A. Unstrained artificial graphene 

It was shown before^ that the two lowest energy bands of 
artificial graphene are similar to the bandstructure of normal 
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FIG. 3: Bandstructure of the unstrained (a = 1) artificial graphene 
for various values of antidot radii p and strengths V of the repulsive 
potential. The bands are shifted by a constant energy to set the energy 
at the Dirac point equal to zero. 


graphene. They touch each other at two Dirac points located 
at the corners K of the hexagonal Brillouin zone. In the vicin¬ 
ity of the Dirac points, the bands have linear dispersion and 
form Dirac cones. Our numerical calculations are in accor¬ 
dance with these results as shown in Fig.[3]for a wide range 
of potential parameters. The four graphs correspond to vari¬ 
ous values of antidot radii p and each graph contains data for 
three values of potential strength V. The free-electron band- 
structure (dashed lines) is shown for comparison. Let’s focus 
on the two lowest energy bands in the M-K direction. One can 
observe that when applying the antidot potential, the two low¬ 
est energy bands (which are doubly-degenerated for V = 0) 
repel each other except at the K point where the Dirac cone 
is formed. This process was analytically described in Ref. [s] 
for small potential strengths in the nearly-free electron (NFE) 
approximation and its generalization for the strained artificial 
graphene is presented in the following section. 

Figure [4j a) shows the ratio of the numerically calculated 
Fermi velocity vf (the slope of the bands crossing at K) and 
the NFE prediction^ vq = 2i:h/ {3m* L) as a function of the 
potential strength for various antidot radii. 

As expected, the NFE approximation is quite accurate for 
small values of the potential strength independent of p/L. 
With increasing potential strength, the Fermi velocity de¬ 
creases, an effect that is more pronounced for larger antidots. 
In the limit V/Eq —> oo, the Fermi velocity saturates at a 
constant value which is relatively close to vo, for p/L < 0.3. 

In order to have Dirac cones as the sole energy values in the 
vicinity of the Dirac points, the energy gap at the M point (see 
Fig .0 should be large enough. This energy gap increases with 
increasing potential strength and for V > V m m, the value of 
the second band at the point M is higher than the Dirac point. 
Then, the Dirac cones are the only energy values around Dirac 
points. Fig. 0|b) shows the minimal value of the potential 
strength V m i n as a function of the antidot radius. For p/L > 



FIG. 4: (a) Ratio of the numerically calculated Fermi velocity vf and 
the NFE value vo as a function of the potential strength for various 
antidot radii p/L. (b) Minimum value of potential strength V m i n to 
obtain only Dirac cone energies around the Dirac point as a function 
of antidot radius p/L. 


0.1, the minimal value of the potential is of the order of Eq 
but for very small p/L the potential must be much stronger. 


B. Artificial graphene under uniaxial strain 

As was shown in section II, the tight-binding model on 
a hexagonal lattice displays Dirac points at the corners of 
the Brillouin zone Ki and Kj. When we apply a hopping 
anisotropy in the armchair direction, the Dirac points are 
shifted and lie somewhere at the line CC / (see Fig. [TJb)) de¬ 
pending on the value of anisotropy parameter t'/t. If t'/t < 1 
the Dirac points are located inside the Brillouin zone and for 
1 < t'/t < 2 the Dirac points lie at the edge of the Brillouin 
zone. For t'/t = 2 both Dirac points merge at the point S 
and for t'/t > 2, an energy gap appears. The aim of this sub¬ 
section is to show a similar behavior in the case of electronic 
artificial graphene. 

In normal graphene, the hopping anisotropy is induced by 
applying a mechanical strain. When the lattice is stretched 
(or compressed), then t'/t < 1 (or t'/t > 1). However, to 
describe a realistic model of normal graphene, the hexagonal 
model from section II is not sufficient because the change of 
the lattice geometry due to the applied strain must be taken 
into accountJ^ial But the results of the realistic modeU- cap¬ 
ture essentially the same features as in the case of the unde¬ 
formed hexagonal model discussed above. 

In the case of electronic artificial graphene, the effect of 
strain can be obtained by tuning the lattice constant in the 
armchair direction. According to ©. the values 0 < a < 1 
correspond to stretching and 1 < a correspond to compres¬ 
sion. Like in normal graphene, the lattice deformations imply 
the deformations of the Brillouin zone. When 0 < a < v/3, 
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FIG. 5: Bandstracture (solid lines) of artificial graphene changing 
with the applied uniaxial strain a. The bands correspond to the CC' 
line in the repeated Brillouin zone scheme (see Fig. [2}b)). Figures 
(a) and (b) show the case of strained artificial graphene, figure (c) 
corresponds to the unstrained case and figures (d)-(f) depict the case 
of compressed artificial graphene. The dashed lines show the posi¬ 
tions of the corners of Brillouin zone Ki and K( which change with 
applied strain according to 1121 . The results are calculated for the 
potential strength V/Eo = 10 and antidot radius p/L = 0.3. 


the Brillouin zone is hexagonal, but the length of its sides to¬ 
gether with the position of the K points depend on a as shown 
in (fTTb and (THT i. 

Figure 0 shows the bandstructure (solid lines) of artificial 
graphene under uniaxial strain for potential strength V/Eq = 
10 and antidot radius p/L = 0.3. For simplicity, we show 
only the bands on the CC / line which captures the positions 
of both Dirac points (see Fig. EJb)). The graphs correspond 
to various values of strain a and Fig. d c ) represents the un¬ 
strained case where a = 1. Here, both Dirac points lie at 
the Ki and K’ i points depicted by dashed lines. In the case 
of stretching (Fig. [5ja), (b)), the Dirac points are located in¬ 
side the Brillouin zone (outside the K i K\ edge) and with in¬ 
creasing stretching (decreasing a) the Dirac points are shifted 
towards the ends of CC ! line. For a = 0.5 the Dirac points 
almost coincide with C and C / . On the other hand, in the case 
of compressive strain where a > 1 (Fig. HJd)-(f)), the Dirac 
points lie at the edge K | K/ and with increasing compression 
(increasing a ) they approach each other. For a = 1.17 the 
Dirac points merge at the point S and the dispersion becomes 
parabolic as shown in Ref. [lfl For a > 1.17, the system un¬ 
dergoes a topological transition, Dirac points no longer exist 
and the energy gap appears (Fig. [5jf)). 

To study the positions of the Dirac points as a function of 
strain, we calculate the distance 7 of the Dirac points from the 
S point given by the formula 


|Dt ~ S| |P 2 -S| 
|C-S| |C' — S| ’ 


(19) 



a a 

FIG. 6: Relative position of the Dirac points with respect to the S- 
point (see Fig. |2}b)) as a function of strain. The graphs correspond 
to the various antidot radii and contain the data for different potential 
strengths and these parameters reflect those from Fig. [3] The dashed 
lines depict the position of the K points given by 01. 


where |C — S| = |C' — S| = 7 t/L. If 7 = 1, the Dirac points 
lie at C and C' and if 7 = 0 the Dirac points meet and merge 
at the S point. 

Figured shows the positions of the Dirac points as a func¬ 
tion of strain for four various antidot sizes and three poten¬ 
tial strengths corresponding to the parameters in Fig. d The 
dashed lines show the position of the points Ki and K( with 
respect to S given by (IT2l) as 


|Ki — S| |k; - S| a 2 

|C — S| |C'-S| 3 


( 20 ) 


For a = \/3 ~ 1.73, the points Ki, K/, S merge ( 7 k = 0) 
and the antidot lattice becomes square instead of triangular. 

Figure d confirms that the behavior described in Fig. d is 
universal: In the case of stretched artificial graphene (a < 1), 
the Dirac points are located inside the Brillouin zone with the 
asymptotic position at C and C' for a —> 0. In the case of zero 
strain (a = 1) the Dirac points lie at Ki and K( with 7 = 2/3. 
When the strain is compressive (a > 1), the Dirac points are 
located at the edge K 1 K/ and move towards each other. At a 
critical value of strain parameter a c , the Dirac points merge 
at S. For a > ac the energy gap increases with increasing 
strain. 

Figure d shows how ac depends on the antidot radius for 
various values of potential strength. For large and strong po¬ 
tential, ac only slightly differs from 1 which means that the 
artificial graphene undergoes the topological transition when 
applying very small compressive strain. On the other hand, 
in the limit of small antidots, the critical strain parameter is 
close to y/3 for arbitrary potential strength. This means that 
the Dirac points are always located at the corners of the Bril- 
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FIG. 7: Critical strain ac corresponding to the merging of the Dirac 
points as a function of the antidot radius for various values of poten¬ 
tial strength V. 


is a dimensionless antidot volume. We can replace the term 
Ji{Gjp)/{Gjp) in (l22l) by its maximum value 0.5 (for p/L —> 
0) and finally obtain the condition for the validity of the NFE 
approximation in the form 


a -C 


3^3 

27T 


1. 


(24) 


A. Dirac points 


Now we are about to derive the formula for the position 7 
of the Dirac points as a function of strain. 

According to the NFE approach^ we can express the elec¬ 
tron wavefunction for k in the vicinity of the Dirac point D 
(see Fig. Ob)) in the basis 


piDi r 

111 = |2> 


as 


i(Di-aJ)-r i(Di-aJ)-r 

7£ ’ |3)= 7£ ’ 

(25) 


louin zone and the system undergoes the topological transition 
from a conductor into an insulator together with the transition 
from triangular to square lattice. 


V. ANALYTICAL MODEL FOR WEAK POTENTIALS 


In this section our goal is to explain analytically the results 
from the previous chapter on the basis of the NFE approach 
developed in Ref. 0 To calculate a first order correction to 
the free-electron energies, the potential has to be sufficiently 
small. Therefore first we examine the validity of the NFE 
approach with respect to the potential parameters. 

For simplicity, we assume the unstrained artificial 
graphene, where the Dirac points are located at K-points and 
the Brillouin zone has a shape of a regular hexagon (the same 
as in Fig. [TJb)). We are interested in the lowest energy band 
around the Dirac point located at the point Ki. In the NFE 
approach, the lowest-energy electron wavefunction (IT5l > for 
k close to Ki can be approximated as a sum of three plane 
waves e 4 k_G j )'7 where G, = {0, a). a 2 } with free-electron 
energies close to Eq given by (IT8l) . The next most significant 
plane waves in the expansion (IT5l) would have energies close 
to 4F 0 . Therefore, the NFE approach is valid when 


\V Gj \C4E 0 -E 0 , (21) 


VK r) = e iqr 
where we have used k = 

Di = 


[6 1 |1)+6 2 |2)+6 3 |3)] i 
q + Di and 


(26) 


(27) 


is the position vector of the Dirac point Di. Then we can write 
the Hamiltonian ( [T~4l > in the basis ( |251 ) as 


H = H 0 + H 1 , 


(28) 


where 


Hn = 


( 2FH D l| 

V 


and for qL <C 1 


14t 


Kf-al 


*|2 


h* 
2m* 


14;-s 

|Di — a 2 


*12 


(29) 



2q-Di 0 0 \ 

0 2q • (D x — a*) 0 

0 0 2q (Di a 2 ) J 

(30) 


In d29l >. we have omitted the diagonal terms 14 which only 
cause a shift of the energy eigenvalues. Because |a*| = |a 2 | 
and | D 1 — aj | = | D 1 — a 2 1 the Hamiltonian Ilf, can be written 
in the form 


where Vg is given by (fTTI) . Condition d2TT> can be rearranged 
to 


47T 1 

v/3 Gjp 


J, (Gjp) <C 3, 


where 




2 


( 22 ) 


(23) 


( Wi Ai A 1 \ 

Ho = [ A\ W 2 A-2 I , (31) 

y A\ A2 W2 J 


where Ai = V a * = V - a * = V a * = V- a *. A 2 = 14;;-a* = 
K*-a;, W\ = ^IDil 2 and W 2 = - a ^| 2 = 

We are searching for the position of the Dirac point D\ and 
therefore we expect that for q = 0 two energy bands touch at 
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D i. This means that the Hamiltonian H (q = 0) = Ho should 
have a pair of equal eigenvalues. The eigenvalues of ( |3T1 > are 

Ei = ~A 2 + IT 2 ; 

£2 = ~ ^^2 + Wi + W 2 — + (A 2 4- W 2 — Wi) 2 ^ , 

63 = — ^^-2 + Wi + W 2 + \J 8 A 2 + (A 2 4- W 2 — Wi) 2 ^ . 

We assume that the first and second eigenvalue are equal 
which yields the condition 

A 2 -A 2 + A 2 {W 2 - Wi) = 0. (32) 


This condition determines the position 7 of the Dirac point Di 
which is hidden in the diagonal terms W\ and W 2 in the form 


W 2 - Wi 



(33) 



We first explore the limit p/L -C 1. In this case the Fourier 
coefficients V p given by O can be simplified as 


Vq = 


27rV p 
~S~Q 


Ji{Qp) 


27tV p (Qp\ 

~S~Q \Y ) 


7rp 2 V 

S 


= V 0 , 
(34) 


which is independent on the wavevector Q. Therefore Ai = 
A 2 and from d32l ) we get W\ = W 2 . Then, from ( l33l > we 
get 7 = 1 — a 2 /3 which is exactly the same as (f20t . This 
means that for p/L < 1 the Dirac point D\ is located di¬ 
rectly at the corner K\ of the Brillouin zone. This explains 
the behavior observed in Fig. [ 6 j where with decreasing anti¬ 
dot radius the position of Dirac points is closer to the dashed 
line given by (l20l >. For p/L = 0.1 and V/Eq = 2, the posi¬ 
tion of Dirac points and K-points is indistinguishable. Then 
the critical strain ac, which leads to the merging of the Dirac 
points, is equal to the strain corresponding to the transition 
from triangular to square antidot lattice, i.e., ac = x/3- 
Now we put the expressions (fl7l ) and ( |3H into (l32t and 
obtain 7 as a function of strain a and relative barrier size p/L 
in the form 



(35) 


where cr is the dimensionless volume of the barrier given by 

<E3. 

Fig.[ 8 ]shows the position of the Dirac points 7 as a function 
of the strain parameter a for weak potentials with the antidot 
volume a = 0.08. The numerical data (filled symbols) and 
corresponding analytical formula (l35l > (solid lines) are plotted 
for various antidot radii shown in the legend. The analytical 
formula roughly agrees with the numerical data except for the 
values of a in the vicinity of critical value ac- Here, the NFE 
basis ( |25| > is not sufficient because there exists another plane 
wave with similar energy (for G = aj + a;)). 


FIG. 8: Position of the Dirac points with respect to the point S as 
a function of the strain parameter for weak potentials with constant 
antidot volume a = 0.08 and various radii listed in the legend. The 
numerical data are shown as filled symbols and solid lines correspond 
to the analytical formula 1351 . The dashed line represents the position 
of the K-point given by @0). 


VI. ABOUT THE EXPERIMENTAL REALIZABILITY 

In experiments with artificial graphene*^ the merging of 
the Dirac points is usually revealed when measuring the den¬ 
sity of states (or some related quantity). When the Dirac 
points exist, the DOS is gapless and equal zero at the energy 
of the Dirac point. When the Dirac points merge, the energy 
gap opens which is visible also in the DOS. In this section, we 
propose the experimental observation of Dirac point merging 
in the electronic systems based on the artificial graphene stud¬ 
ied in this paper. 

We first discuss the conditions under which the energy gap 
in the DOS opens as a result of Dirac point merging. To ob¬ 
serve the energy gap in the DOS, no electron energies should 
lie in the energy window given by the gap between the bands. 
Otherwise, the energy gap in the DOS would not open even if 
the Dirac points merged. To analyze this problem, it is suffi¬ 
cient to study the bandstructure at the edges of the Brillouin 
zone, because here the maximum of the first band and mini¬ 
mum of the second band are located. Fig. [9] shows how these 
bands change with increasing compressive strain. The graphs 
represent the energy bands corresponding to the MKiS line at 
the edge of the Brillouin zone (see Fig.[2jb)) calculated for an¬ 
tidot radius p/L = 0.2 and two values of the potential strength 
V/E 0 = 3 (solid lines) and V/Eq = 6 (dashed lines). 

Fig- 0a) shows the data for unstrained artificial graphene. 
According to Fig. 0B), both potential strengths satisfy the 
condition V > V m i n which means that a pair of Dirac cones 
located at Ki and K ’1 (not shown here) are the sole energy 
values in the vicinity of Dirac point. With increasing com¬ 
pression (Figs. Hfb)-(d)), one can observe at the right part of 
each panel that the Dirac point moves to the point S where it 
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FIG. 9: The development of the bandstructure of artificial graphene 
with increasing compression given by the strain parameter a shown 
in the legend. The bands correspond to the line MKiS in the Bril- 
louin zone (see Fig.|2jbj) where the maximum of the first band and 
the minimum of the second band are located. The parameters are 
p/L = 0.2 and V/Eo = 3 (solid lines) or V/Eo = 6 (dashed 
lines). Graph (a) shows the bands of unstrained artificial graphene 
where the Dirac cones are the only energy values around the Dirac 
points. Graphs (b) and (c) depict the merging of the Dirac point for 
V/E 0 = 6 and V/Eo = 3, respectively. The graph (d) shows the 
situation where the Dirac points are already merged and the gap be¬ 
tween the bands exists. In the case of V/Eo = 6, the maximum of 
the first band is lower than the minimum of the second band, which 
means that the gap appears also in the DOS. In the case of V/Eo = 3, 
the gap between the bands is obscured by the second energy band and 
the DOS remains gapless. 


merges with the second Dirac point (not shown) and the gap 
opens. At the same time, the energies of the second band in 
the left part of the panels are getting lower with respect to 
the Dirac point. When the Dirac points merge for V/Eo = 3 
(Fig- HI c)), the energy minimum of the second band located at 
the M point is already lower than the Dirac point. Therefore, 
when the gap in S opens (Fig. |9jd)), it does not appear in the 
the DOS because it is obscured by the second energy band. 
On the other hand, in the case of stronger potential V/Eq = 6 
when the Dirac points merge (Fig. |9jb)), the energies of the 
second band are well above the Dirac point. When the energy 
gap opens (Figs. |9jc,d)), the energies of the second band are 
still higher than the maximum of the first band located at the 
point S. This means that the energy gap is visible also in the 
DOS. 

The above analysis of the bandstructure has revealed that 
for every p/L there is a minimal value of potential strength 
V 1 [ 1 j ] ( t ) 1 ,s that is necessary to observe the gap opening in the DOS 
and this value is larger than V m i n given by Fig.[4] b). We obtain 
Fin in S numerically using the constraint that for critical strain 
q/DOS w j len th e Dirac points merge, the energy value of the 



FIG. 10: (a) The minimal value of the potential strength that 

is needed to observe the gap in the DOS. This value is calculated as 
a function of the normalized antidot radius p/L. To obtain V°£ S , 
two conditions have to be satisfied simultaneously. Firstly, the Dirac 
points merge and the bands touch each other only at the point S. 
Secondly, the minimum of the second energy band, which is located 
at the point M (see Fig. m, should be equal to the energy of the Dirac 
point. This means that for V = V^,° s the minimum of the second 
energy band is equal to the maximum of the first band, (b) Critical 
strain corresponding to the potential with the strength V^,° s . 


second band at the point M should be equal to the energy of 
the Dirac point. Then, for the potentials larger than 
there is an interval of strain parameter a where the gap in the 
DOS appears. For V < the gapped DOS cannot be 

observed regardless of the value of strain a. Fig. m shows 
the calculated together with the corresponding critical 

strain a° os . 

There have been a couple of attempts to observe the Dirac 
quasiparticles in the experimental electronic systems based on 
the theoretical model of artificial graphene studied in this pa¬ 
per. In the work of Ref. [ToL the authors discussed a possible 
realization of artificial graphene by the two-dimensional elec¬ 
tron gas (2DEG) in a semiconductor heterostructure where 
the triangular lattice of antidots is created by etching or gat¬ 
ing. Moreover, they performed an experiment with a 2DEG 
made from GaAs (with the effective mass of the electrons 
to* = 0.067m e ) where the antidot lattice was etched with the 
lattice constant L = 200 nm and p/L ~ 0.15. The estimated 
potential strength was about 2-4 meV. These parameters cor¬ 
respond to V/Eo ~ 10 where Eq ~ 0.25 meV. Although the 
presence of Dirac electrons was not fully confirmed in their 
system (due to the strong disorder) the authors conclude that 
this should be technologically feasible by reducing the lattice 
constant of the antidot lattice below 100 nm. 

The same concept of artificial graphene was used in the 
experiments with so-called molecular graphene^ 7 where the 
2DEG on the copper surface is modulated by the triangular 



























10 



E [meV] E [meV] 



E [meV] E [meV] 


FIG. 11: The DOS of electrons in artificial graphene calculated as 
a function of the energy for various values of strain parameter a 
shown in the legend. The parameters p = 0.7 nm, L = 2.7 nm, 
and V = 1 eV used in the simulation are the same as in the experi¬ 
ment with molecular graphene created by coronene! 7 (a) corresponds 
to the strained artificial graphene, (b) shows the unstrained case, and 
(c)-(d) corresponds to the artificial graphene compressed in the arm¬ 
chair direction. The gray region in (d) depicts the energy gap created 
as a result of the merging of the Dirac points. 


antidot lattice of repulsive molecules (CO or coronene). The 
manipulation of several hundreds of molecules is provided by 
a scanning tunneling microscopy. The effective mass of sur¬ 
face electrons is m* = 0.38m e and the lattice constant of 
the antidot lattice L = 1 — 3 nm varied for different exper¬ 
imental situations. The radius of the antidot depends on the 
details of the electrostatic potential of the repulsive molecule 
and in the case of coronene it is estimated^ to be p = 0.7 nm. 
In the case of CO, the authors of Ref. |h| do not provide the 
exact value of p but it can be roughly estimated from the to¬ 
pographs depicting the distribution of electron density in the 
modulated 2DEG as p ~ 0.4 nm. A potential strength of 
an antidot was estimated roughly to be V ~ 1 eV for both 
types of molecules. The presence of Dirac quasiparticles was 
demonstrated by the measurements of the differential conduc¬ 
tance between the STM tip and the Cu surface that mimic the 
DOS of modulated electrons. The spatially averaged differen¬ 
tial conductance spectrum showed the features typical for the 
DOS of graphene: a zero value at the Dirac point, linear shape 
around it and a pair of van Hove singularities. 

For the purposes of this paper, particularly interesting is 
the experiment described in the supporting online material 
of the work in Ref. [b| where the authors studied the effect 
of uniaxial strain on CO molecular graphene. Here, the lat¬ 
tice constant was L = 1.785 nm which gives V/Eq ~ 2 
and p/L ~ 0.4. The measurement of the differential conduc¬ 
tance spectrum showed that when the molecular graphene is 
elongated (stretched) by about 30% in the armchair direction 
the DOS around the Dirac point remains linear and gapless. 


This is in agreement with our calculations which predict that 
in the stretched artificial graphene the Dirac points exist and 
they are shifted inside the Brillouin zone. We would like to 
stress that according to our results, the merging of the Dirac 
points appears if the artificial graphene is compressed (instead 
of stretched) in the armchair direction. 

According to Fig.|T0ja), the potential strength V/Eq ~ 2 
at p/L ~ 0.4 is close to the value V°P S . Therefore, it is not 
sure that it would be possible to observe the merging of Dirac 
points in the experiment with CO molecules. More promis¬ 
ing is the molecular graphene made from coronene.? 7 where 
the lattice constant is L = 2.7 nm, which gives p/L = 0.26 
and V/Eq = 4, a value well above V°° s . Fig.Qj]shows the 
DOS of artificial graphene calculated for these experimental 
parameters. The different graphs show how the DOS changes 
with strain. Figure fTTTb) shows the DOS for unstrained arti¬ 
ficial graphene, which is very similar to the DOS of normal 
graphene. The graph|TT]a) corresponds to the situation where 
the artificial graphene is strained by about 30% which mod¬ 
els the experimental situation discussed above. As expected, 
the DOS is gapless because the Dirac points do not merge. On 
the other hand, the granhs UTT c-d) show how the DOS changes 
with compression. The merging of the Dirac points occurs for 
ac = 1.33, therefore for a = 1.5 (Fig.QlJd)), the gap in the 
DOS (gray region) is nicely recognizable. The value a = 1.5 
corresponds to the compression of the antidot lattice of about 
100(1 — 1/a) ~ 33% which should be possible to realize in 
the experiment with molecular graphene. 

As discussed in this paper, merging of the Dirac points oc¬ 
curs in artificial graphene under uniaxial compression in the 
armchair direction. We would like to stress that the same ef¬ 
fect can be obtained by stretching the antidot lattice in the 
zigzag (horizontal) direction. In our model, stretching in the 
zigzag direction just increases the lattice constant L to L" 
while the vertical dimension of the antidot lattice remains the 
same. Therefore, the antidot lattice stretched in zigzag direc¬ 
tion is the same as the lattice with lattice constant L" com¬ 
pressed in the armchair direction. 


VII. CONCLUSIONS 

In normal graphene, the theory based on tight-binding ap¬ 
proximation predicts that with increasing anisotropy in the 
hopping matrix elements, both Dirac points are moving along 
the edge of the Brillouin zone towards each other until they 
merge. This anisotropy can be realized by application of com¬ 
pressive strain in the armchair direction. A merging of the 
Dirac points was observed so far only in artificial systems that 
mimic the properties of graphene: in the experiments with 
confined microwaves in a hexagonal array of waveguides^ 
and in a laser lattice with cold atoms?2J2 

We have shown numerically and analytically that the merg¬ 
ing of the Dirac points can be observed also in electronic arti¬ 
ficial graphene. The artificial graphene we considered is to be 
created from the two-dimensional electron gas by applying a 
repulsive triangular potential and the effect of strain was mod¬ 
eled by reducing the distance between the repulsive potentials 
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along the armchair direction. 

Our numerical calculations have show that molecular 
graphene with coronene^ is a promising candidate to observe 


the merging of the Dirac points. This should occur for a lat¬ 
tice compression of about 25%, which appears technologi¬ 
cally feasible. 


* Electronic address: juraj.feilhauer@ptb.de 

1 A. H. Castro Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and 
A. K. Geim, Rev. Mod. Phys. 81, 109 (2009). 

2 M. Polini, F. Guinea, M. Fewenstein, H. C. Manoharan, and V. 
Pellegrini, Nat. Nanotechnol. 8, 625 (2013). 

1 L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, 
Nature 483. 302 (2012). 

4 M. Bellec, U. Kuhl, G. Montambaux, and F. Mortessagne, Phys. 
Rev. Lett. 110, 033902 (2013). 

5 Y. Plotnik et al., Nat. Mater. 13, 57 (2013) 

6 K. K. Gomes, W. Mar, W. Ko, F. Guinea and H. C. Manoharan, 
Nature 483, 306 (2012) and corresponding supplemental online 
material. 

7 S. Wang, L. Z. Tan, W. Wang, S. G. Fouie and N. Fin, Phys. Rev. 
Lett. 113, 196803 (2014). 

8 C.-H. Park and S. G. Fouie, Nano Lett. 9, 1793 (2009). 

9 M. Gibertini, A. Singha, V. Pellegrini, M. Polini, G. Vignale, 
A. Pinczuk, L. N. Pfeiffer, and K. W. West, Phys. Rev. B 79, 
241406(R) (2009). 

10 L. Nadvormk, M. Orlita, N. A. Goncharuk, F. Smrcka, V. Novak, 
V. Jurka, K. Hruska, Z. Vyborny, Z. R. Wasilewski, M. Potemski, 
and K. Vyborny, New J. Phys. 14, 053002 (2012). 

11 E. Kalesaki, C. Delerue, C. Morais Smith, W. Beugeling, G. Al¬ 
lan, and D. Vanmaekelbergh, Phys. Rev. X 4, 011010 (2014). 

12 Y. Hasegawa, R. Konno, H. Nakano, and M. Kohmoto, Phys. Rev. 


B 74, 033413 (2006). 

13 V. M. Pereira, A. H. Castro Neto, and N. M. R. Peres, Phys. Rev. 
B 80, 045401 (2009). 

14 P. Dietl, F. Piechon, and G. Montambaux, Phys. Rev. Lett. 100, 
236405 (2008). 

15 J.-N. Fuchs, arXiv: 1306.0380 

16 G. Montambaux, F. Piechon, J.-N. Fuchs, and M. O. Goerbig, 
Phys. Rev B 80, 153412 (2009). 

17 S.-L. Zhu, B. Wang, and L.-M. Duan, Phys. Rev Lett. 98, 260402 
(2007). 

18 B. Wunsch, F. Guinea, and F. Sols, New J. Phys. 10, 103027 
(2008). 

19 O. A. Tkachenko and V. A. Tkachenko, JETP Lett. 99, 204 (2014). 

20 O. P. Sushkov and A. H. Castro Neto, Phys. Rev. Lett. 110, 186601 
(2013). 

21 C.-H. Park, L. Yang, Y.-W. Son, M. L. Cohen, and S. G. Louie, 
Phys. Rev Lett. 101, 126804 (2008). 

22 M. Aichinger, S. Janecek, I. Kylanpaa, and E. Rasanen, Phys. Rev 
B 89, 235433 (2014). 

23 Z. Liu, J. Wang, and J. Li, Phys. Chem. Chem. Phys. 15, 18855 
(2013). 

24 M. R. Masir, D. Moldovan, and F. M. Peeters, Solid State Com- 
mun. 175-176, 76 (2013). 



